---
title: "Roads"
author: "Francesco Bailo"
date: "23/07/2020"
output: html_document
---

```{r eval=F}
library(sf)

road_classes <- 
  c("motorway", "motorway_link",
    "primary", "primary_link",
    "trunk", "trunk_link")
  
ita_osm_roads.sf <- 
  rbind(sf::st_read("/pvol/gis/ita_osm_nordest/gis_osm_roads_free_1.shp") %>%
          dplyr::filter(fclass %in% road_classes),
        sf::st_read("/pvol/gis/ita_osm_nordovest/gis_osm_roads_free_1.shp") %>%
          dplyr::filter(fclass %in% road_classes),
        sf::st_read("/pvol/gis/ita_osm_centro/gis_osm_roads_free_1.shp") %>%
          dplyr::filter(fclass %in% road_classes),
        sf::st_read("/pvol/gis/ita_osm_sud/gis_osm_roads_free_1.shp") %>%
          dplyr::filter(fclass %in% road_classes),
        sf::st_read("/pvol/gis/ita_osm_isole/gis_osm_roads_free_1.shp") %>%
          dplyr::filter(fclass %in% road_classes))

save(ita_osm_roads.sf, file = "../roads/ita_osm_roads.sf.RData")
```


```{r}
library(sf)
library(parallel)
library(dplyr)

load("../roads/ita_osm_roads.sf.RData")
load("../grid/hex_ita_10km.sf.RData")

ita_osm_roads.sf <- 
  st_transform(ita_osm_roads.sf, 32632)

hexRoadTotLength <- function(i, hex_ita_10km.sf, ita_osm_roads.sf) {
  require(sf)
  this_hex <- 
    hex_ita_10km.sf[i,]
  these_roads <- 
    sf::st_intersection(ita_osm_roads.sf, this_hex)
  this_length <- 
    sum(as.numeric(sf::st_length(these_roads)))
  return(data.frame(hex_id = this_hex$hex_id,
                    tot_length = this_length))
}

cl <- makeCluster(10)
res.list <- 
  parLapply(cl, 1:nrow(hex_ita_10km.sf), hexRoadTotLength, 
            hex_ita_10km.sf, ita_osm_roads.sf)
primary_road_hex.df <- bind_rows(res.list)
stopCluster(cl)

save(primary_road_hex.df, file = "../roads/primary_road_hex.df.RData")

```


```{r}
load("../roads/primary_road_hex.df.RData")

hex_ita_10km_roads.sf <- hex_ita_10km.sf

ggplot

hex_ita_10km_roads.sf$road_length <- 
  primary_road_hex.df$tot_length[match(hex_ita_10km_roads.sf$hex_id,
                                       primary_road_hex.df$hex_id)]


ggsave(filename = "../fig/osm_primary_road.png", width = 12, height = 9,
       plot_grid(
         ggplot(ita_osm_roads.sf) +
           geom_sf(size = .1, colour = 'black') +
           theme_bw() +
           theme(axis.text = element_blank(),
                 axis.ticks = element_blank()) +
           labs(title = "OpenStreeMap primary road network"),
         ggplot(hex_ita_10km_roads.sf %>% st_transform(4326)) +
           geom_sf(aes(fill = road_length), colour = NA) +
           scale_fill_distiller(labels = function(x)paste0(format(x, big.mark=",")," km"), 
                                palette = "Reds", direction = 1) +
           theme_bw() +
           theme(axis.text = element_blank(),
                 axis.ticks = element_blank()) +
           labs(fill = NULL, title = "Total road length by hexagon cell"),
         ncol = 2, align = "h"))

ggplot(hex_ita_10km_roads.sf) +
  geom_density(aes(road_length))

```

